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Abstract 


In this wurk. an axisyiinm'tric large' dofonnatioii elasto-plastic finite element code for deter- 
mination of damage in cold forging process is developed. For typical values of input variables, 
th(' damage, hydrostatic stress, circumferential stress and axial stress distributions are stud- 
i('il. A detailed parametric study of damage is carried out to show the effects of the coefficient 
of friction at the die-work interface and the height to diameter ratio of the blank. 

The Updated Lagrangian formulation, which is convenient for handling material and ge- 
ometric non-linearitie.s, is used. New incremental objective stress measure and logarithmic 
strain measure are used which allow the use of large incremental displacement. The forging 
process is a displacement-control problem. Therefore, to accelerate the convergence of the 
iterative solution scheme, arc length method is used in conjunction with the modified Newtot 
Raphson iterative technicpie. The material is assumed to be elastic-plastic strain hardening 
and yit'lding accmrdiiig to the von Mises criterion. Coulomb’s friction law is used to mode 
the friction at the die-workpiece interface. To facilitate the application of friction boundarj 
condition and also to reduce the computational time, the static condensation of the coeffi> 
cient matrix is carried out. The inertial and body forces are neglected. A continuum damagf 
mechanics criterion to model ductile fracture in metal forming is used for defect prediction. 

The study shows that the micro-cracks first initiate along the die-workpiece interfact 
region near to the free surface, then at the centre of the workpiece and proceed to the meridiaj 
surface. The fracture will appear first at the meridian surface since, here, the hydrostatic stresj 
and axial stress are the minimum compressive and the circumferential stress is the maximuri 
tensile. Thus, micro-cracks here can grow faster to a fracture. From the parametric studj 
it is observed that at higher friction, micro-cracks initiation and meridian surface fracturi 
occurrence take place at less reduction. Bulging is more severe when the friction at die 
work interface is higher. As compared with the case of height to diameter ratio equal t| 
one, blanks with height to diameter ratio less than or greater than unity can sustain moil 
reductions without micro-cracks and fracture. Further, with higher height to diameter rati(| 
the chances of formation of a central cavity are more because of less compressive hydrostat| 
stress at the centre. I 
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Chapter 1 


INTRODUCTION AND 
LITERATURE REVIEW 


1.1 Introduction 

Forging is a process in which the workpiece is shaped by compressive forces applied through 
various dies and tools. It is one of the oldest metalworking processes known to the mankind. 
It was first used to make jewelry, coins and various implements. Today, from a simple rivet 
to the landing gear of aircraft or from an automobile connecting rod to the propeller shaft 
for ship, forged products have quite a wide spectrum of applications. Good strength and 
toughness of the forged parts is the main reason for this. These properties can be achieved 
since grain orientation and grain size are controllable in forging process. 

Forging may be performed at room temperature (cold forging) or at elevated tempera- 
tures (hot forging). Cold-forged parts have good surface finish and dimensional accuracy 
but because of higher strength of the material, cold forging requires greater forces, hence 
the workpiece materials must have sufficient ductility at room temperature. Less ductility, 
micro-structural faults in original metal, incorrect die design, improper heating and lubrica- 
tion, large and rapid reduction, nonuniform deformation of material etc. can cause forging 
defects. Barreling, tears (surface cracks), edge cracks, central cavity, thermal cracks, laps, 
fins, wrinkles (orange peel), end grains etc. are the defects seen in forged parts [54]. Defects 
are undesirable because these reduce the strength of the forged parts which are intended to be 
used in critical stress applications. Several ductile fracture criteria are used to predict cracks 
in forging [65,69] and need exact stress and strain values to predict micro-crack formation. 
The finite element method can estimate stresses and strains to a good degree of accuracy. In 
the present work, an elasto-plastic FEM analysis of axisymmetric upsetting of cylindrical bil- 
let is carried out and damage distribution is predicted using a continuum damage mechanics 
model, proposed recently by Dhar [46,75]. 



1.2 Literature Review 

1.2.1 Deformation Analysis of Forging Process 

The problem of upsetting of metals has theoretical as well as technological significance in 
that its solution help in bringing about a better understanding of the behaviour of metals 
during pla.stie deformation at high pressure. Such solutions may also give specific answers to 
problems actually encountered in practice. It is for this reason that a number of investigators 
have dealt with the forging problem. Early investigators considered the material as an ideal- 
plastic solid i.e. rigid perfectly-plastic solid. This leads to so called slip line solution. Such 
solutions have been proposed by Prandtl, Hill et al. [1], and Green [2]. The slip line theory 
has been well d(weloped to analyze a non-homogeneous plane strain deformation of an ideal- 
plastic solid. This theory can be used to get very good first approximations to loads required 
or to provide indications of the manner in which material deforms. However, the elastic effect 
can not be incorporated and work hardening is difficult to incorporate in this method and it 
is not suitable for axisymmetric problems. Further, residual stresses also cannot be studied 
using this technique. Shabaik [3] investigated the bulging in upsetting by using the slip line 
theory. 

Hoffman and Sach [4] proposed the slab method. A complete analysis of the slab method 
has been presented by Altan [5] for the axisymmetric closed die forging. In this method, 
a slab of infinitesimal thickness is selected and then the force equilibrium equation for the 
slab is written, assuming that the deformation is homogeneous within the slab and the major 
directions are the principal stress directions. The resulting differential equation is then solved 
with appropriate boundary conditions. While this method can give very good predictions of 
the load variation with deformation, it is inherently incapable of predicting shape changes 
such as barrelling in open die forging. Further, since this method is primarily based on the 
assumption of homogeneous deformation, it possibly cannot predict the residual stresses as 
they are caused by non-uniform deformation. 

The upper bound theorem was formulated by Prager and Hodge [6]. If surfaces of velocity 
discontinuities are included [7], then it states that among all kinematically admissible velocity 
fields V*, the actual one minimizes the following expression 

r = f f T|Au*|ds- [ TiV*ds, (1.1) 

Jv J Sr ^ St 

Here, e*j is the strain rate field derived from v*, S-j is the deviatoric stress field derived from 
e*j and |Au*| is the discontinuity in the tangential component of the velocity across the surface 
Sr- Further, r denotes the shear stress on the surface Sr and T, is the prescribed traction 
on the part St of the boundary. The first term expresses the power spent in causing the 
internal deformation over the volume V of the deforming body. The second term represents 



the power loss over the surfaces of velocity discontinuities including the boundary between 
tool and material. 1 he last term represents the power supplied by specific tractions. What 
the upper bound theorem states is that the actual externally supplied power is never higher 
than that computed by the first two terms of the above equation. 

Kudo [8] applied the upper bound theorem to the problem of plane strain forging and 
extrusion. Kudo [9] also applied this to the analyses of axisymmetric cold forging and ex- 
trusion. McDonald et al. [12], Kobayashi [11], Avitzur [10] and many others have suggested 
upper bound velocity fields to predict forging load in upsetting with bulging. Yang and Kim 
[13], Kim et al. [14] and Manuel et al. [15] have proposed upper bound method for analysis 
of three-dimensional upset forging. The major difficulty in using upper bound method is how 
to choose a kinematically admissible velocity field since the accuracy of the solution depends 
on how close the assumed velocity is to actual one. Another disadvantage of the method is 
that the solution provided by the method does not satisfy the equilibrium equations. So, even 
if elastic effects are included in integral (1.1), the method possibly cannot be used for the 
determination of residual stresses. 

Accurate determination of forging parameters under realistic conditions became possible 
when the finite element method (FEM) was introduced. The major advantage of FEM is that 
it can be applied to a wide class of boundary value problems without restrictions on workpiece 
geometry. Using FEM, it is possible to predict the platen forces, the interface pressure and 
the stress, strain and residual stress distributions within the billet at various deformation 
levels. 

Early applications of the FEM to forging problems were based on the incremental method 
proposed by Lee and Kobayashi [16]. The method uses the elastic-plastic strain matrix 
based on Prandtl-Reuss equations. The additivity of incremental elastic and plastic strains 
is assumed. Even though the stress-strain matrix and the geometry are updated after every 
increment, only the linearized incremental equations are used. The method was applied to 
solid cylinder upsetting [17], ring compression [18] and for predicting defects in upsetting ; 
[19]. Maccarini et al. [20] used the method for studying the influence of die geometry on cold ; 
extrusion forging operation. Where as Lee and Kobayashi [16] used the velocity as the primary | 
unknown. Hartley et al. [21,22] proposed an incremental method with the displacement as | 

the primary unknown. In Hartley’s method, the linearized equations are used and the elasto- J 

plastic matrix and the geometry are updated after every increment. Shima et al. [23] proposed | 
a rigid-plastic finite element model based on plasticity theory for porous materials. They used } 
Coulomb’s law to model the friction at the workpiece-die interface. They applied this method | 
to upsetting of a cylinder and validated their results by conducting experiments. Three- 
dimensional elasto-plastic finite element analysis has been done by Pillinger et al. [24] using 
linearized incremental equations. ! 

Linearized incremeEt4 equatiops give only an approximate solution. If the incremental 



size is not sufficiently small, the error between the exact and the approximate solutions 
grows rapidly with the applied load, as the error in the solution of the current increment 
gets propagated into the next increment when the stress-strain matrix and the geometry are 
updated. To avoid this phenomenon, non-linear incremental equations should be solved by 
using a scheme like Newton-Raphson technique. Such a formulation was first proposed by 
Bathe et al. [25]. this formulation (with or without elastic effects) has been applied to the 
forging problem by quite a few researchers. It was applied by Dadras and Thomas [26] and 
Carter and Lee [27] to axisymmetric upsetting. Dadras and Thomas have also developed the 
kinematically admissible velocity fields for different deformation zones by means of empirical 
equations. 

Some of the typical references in this area are as follows. Most of them use iterative 
method to solve the non-linear incremental equations. Michel and Boyer [28] carried out the 
elasto-visco-plastic finite element analysis of cold upsetting process and did its experimental 
validation. However, they used the friction factor model. Therefore, their results could not be 
used to validate the proposed finite element model, which uses the Coulomb friction model. 

Lin [29] has proposed a new model for inter-facial friction in which the friction coefficient 
changes continuously with the reduction in height. The experimental measurements on the 
barrelling of the free surface are used to determine the functional dependence of the friction 
coefficient on the reduction in height. Further, he assumed the flow stress to depend on strain, 
strain rate and temperature. Even though his model is in good agreement with experiments 
on cold upsetting, it is a little too sophisticated to be of use to practical problems. Yang et al. 

[30] have proposed a 3-D finite element model in which the workpiece and die are analyzed 
together. The workpiece-die contact is classified into three cases: (I) No contact, (II) sliding 
contact and (III) cohesive contact. Multi-substructure technique is used to reduce the size i 
of the problem. Their model incorporates the determination of temperature field as well as 
the residual stress distribution. However, no results on residual stresses are reported. Joun ; 

i 

et al. [31] have proposed a computer simulatioii technique for forging process which uses a | 

spring attached-die. The strategy of spring attached-die control the metal flow lines in such i 

a fashion that it results in prevention of defects and improvement of product quality. ' 

The thrust of most of the above studies has been either to find an accurate estimation } 

of forging load or to study deformation field. A major difficulty observed in all the previous | 

attempts is that the increment size has to be very small in order to achieve a reasonable level | 

of accuracy in the stress values. This leads to prohibitively large amount o computational ■ 

time. The Jauraann stress measure, which is the objective stress measure used presently in ; 

the literature, does not give accurate stress values if the incremental shear deformation is 
large. Thus, there is a need to look for a new objective stress measure, which will allow the 
use of a large incremental size without compromising the accuracy of stresses. Such a measure 
has been recently proposed in reference [33]. 



1.2.2 Ductile Fracture Criteria 


Ductile fracture is often a limiting factor in many metal forming processes. It is well known 
that ductile fracture occurs clue to micro- void nucleation, growth and finally coalescence into 
micro crack. Many ductile fracture criteria have been proposed and applied to predict duc- 
tile fracture in phistically deforming metals including metal forming processes. This section 
reviews some of the published ductile fracture criteria. 

The published ductile fracture criteria can be broadly classified into two groups. They are 

• Models based on microscopic observations about void nucleation, growth and coales- 
cence 

• Empirical and semi-empirical models 

Three broad approaches have emerged which try to predict ductile fracture (micro-crack 
initiation) on the basis of void nucleation, growth and coalescence. They are 

• Models of porous plasticity 

• Models of void nucleation, growth and coalescence 

• Continuum damage mechanics models 

Models of Porous Plasticity 

In these models, the material with voids is idealized as porous material. Thus, its constitutive 
equation is derived from the plastic potential of porous materials. In a porous plastic material, 
there is a possibility of dilatation because of which the yield surface does not remain an 
infinitely long cylinder like that of an incompressible material but is capped by elliptical 
surfaces. 

Gurson’s Model 

Gurson [34] proposed a plastic potential function for a porous solid with randomly distributed 
voids of volume fraction Vf in the form 

*^2 o 

$ = ^ + 2VfCOsh{^) -{1 + V}) = 0 (1.2) 

where dp and amj, are the generalized and hydrostatic stresses of the porous aggregate and a 
is the generalized stress of the void free (incompressible) matrix. 

In this model, the growth rate of void volume fraction is considered partly due to the 
growth of the existing voids and partly due to the nucleation of new voids, i.e., 

Vf = Vf-grouith "b Vf —nucleation (1.3) 



Tho rate of cliangc of the volume fraction due to the growth of existing voids is related 
to the mean strain rate c„ by the relation 

^'f-growth. ~ {1 — Vf)in (1-4) 

In g<’iK’ral. \’t)id iiucleaf ion rale <lopend.s on both the generalized strain rate as well as the rate 
of increase of hydrostatic' stress. In this model, void nucleation rate is assumed to depend 
only on t he generalized strain rate by the relation 

-nucleation ~ (f-^) 

where A is a constant. 

Gurson’s criterion [34] was modified by Tvergaard [35] and Needleman and Tvergaard [36] 
by introducing some constant (gi) to bring the predictions of the model into closer agreement 
with numerical analysis of a periodic array of voids. The modified criterion is given by 

oH 3 (t»« 

* = + ‘2V,q,C03h{^) - (1 + = 0 ( 1 . 6 ) 

where qi is a constant. 

In this model, ductile fracture is regarded as the result of an instability in the dilatational 
plastic flow field localized in a band (called the shear band). The fracture criterion is repre- 
sented as a graph of critical localization strain versus the critical void volume fraction with 
strain hardening exponent as a parameter. 

Oyane’s Model 

Oyane [37] and Oyane et al. [38] used the plasticity theory for porous materials to propose 
a ductile fracture criterion indicating that micro crack initiates whenever the volumetric 
strain reaches a material dependent critical value. For pore free materials, assuming that the 
material obeys porous plasticity equations after the initiation of voids, Oyane obtained the 
criterion for ductile fracture as 

/ = c (1.7) 

^0 

where fp is a function of relative density pr which is defined as the ratio of the apparent 
density of the porous material to the density of the pore free matrix, €„/ is the volumetric 
strain at fracture and n and C are constants. The constant C is given by 

f\l + ^)*- = c ( 1 . 8 ) 

where li is the generalized strain at which voids get initiated, e/ is the generalized strain at 
which fracture occurs and A is a material constant. By assuming that li = 0, the above 
criterion reduces to a simple form 



For porous materials, they modified the criterion by including the relative density term 
into the integral. Thus, 



Here is tlie generalized flow stress of the porous material, S is a constant and po is the 
initial relative density of the porous material. A method for estimating values of the material 
constants using the compression test has been provided by them. 

The major limitation of the above two criterion is that they do not model ductile fracture 
(i.e. micro-crack initiation) as a void coalescence phenomenon. 

Models of Void Nucleation, Growth and Coalescence 

The void nucleation model of Goods and Brown [39], the void growth models of McClintock 
[40] and Rico and Tracey [41] and void coalescence model of Thomason [74] along with some 
additional models on void nucleation and growth are discussed in this section. 

Models of Void Nucleation 

Goods and Browm [39] developed a stress based void nucleation model based on the condi- 
tion that the voids nucleate by decohesion of particles whenever the normal stress reaches a 
critical value ((T^.) at the particle/matrix interface. They obtained a relation between the void 
nucleation strain (e” ) and hydrostatic stress {am) as 

e” = KTj,{ac-amf (1-11) 

where R is a material constant dependent on the particle volume fraction and Vp is the 
particle radius. Experiments on Fe-FesC system [39] confirm the linear relationship between 
the nucleation strain and particle radius. This relation is for small spheroidal particles of 
radius less than or equal to Ipm. 

For large size particles, the approximate analysis of Argon and Im [42] gives the condition 
for the void nucleation by decohesion as 

^ + crm = Oc (1-12) 

For a given problem, this condition can be expressed in terms of the nucleation strain. Note 
that this relation is independent of particle radius. 

Gurland [43] proposed a nucleation model based on the experiments conducted by him on 
1.05% C spheroidal steel. Experimental observations reveal that the void nucleation occur 
by cracking of cementite particles. Voids nucleate at all strain levels depending on the size, 
shape and orientation with the maximum principal stress direction. It is observed that the 
volume fraction of broken particles is a linear function of the plastic strain. 



Void Growth Models 


Once the imcleation of micro-voids takes place either by decohesion or cracking of a secondary 
particle or inclusion, the resulting stress free surface of the void produces a localized stress and 
strain concentration in the plastically surrounding matrix. With continuing plastic flow of the 
matrix, the void undergoes volumetric growth and change of shape. Some of the published 
void growth models are discussed here. 

McClintock [40] proposed a void growth model for two-dimensional plane strain problems 
considering a single elliptic void. He assumed that the major and minor axes of the void 
coincide with the principal stress directions. He obtained the following closed form expression 
for void growth in a rigid work-hardening material whose generalized stress relation is a = 


In 


Rh 


v/Se 


. , \/3 (1 — n) (Ta -H CTj Cfl + €6 

Sinn — + 


(1.13) 


Rqii 2(1 — n) 2 a ' 2 

where Rh and i?o/i are the current and initial mean radii of the hole respectively, Ua and at 
are the principal stress and €a and €5 are the principal strain components along the major and 
minor axes respectively. Fracture was assumed to occur at the point where a growing void 
touches the cell boundary. Dung [44,45] modified and extended McClintock’s model [40] to 
analyze the growth of circular/elliptic voids in plane strain problems and spherical/ellipsoidal 
voids in thrcc-rlimensional problems. 

Rice and Tracey [41] considered a single spherical void of initial radius Rov in a remote 
uniform strain rate field cp and remote stress field aij in a rigid plastic material. They 
obtained the following expression for the rate of change of the radii {Rkv) of the void in the 
principal strain rate directions; 


Rkv = {(1 + T’jejfc -f ean 


(1.14) 


where 

F = I for linear hardening and for low values of am for non hardening 
= 1 for high values of am for non hardening 

H = 0.75^ for linear hardening 

= 0.558 sinh |^ + 0.008z/cosh|^ for non hardening, 

Rmean ~ 3 (-^Iti 4" R2v 4" Rzv) 


i = generalized strain rate 


4 = principal strain rates 

u = -T^^, Lode variable. 


Note that an initial spherical void grows into an ellipsoidal void of principal radii Ri^, R 2 v 
and R^v Thomason [74] integrated the above equation by assuming that the principal axes of 



strain rates remain fixed in direction throughout the strain path and obtained the following 
expressions for the principal radii of the void: 


where 


D - f A t R 

(1.16) 


(1.16) 


(1.17) 

3 + ^ 



and e? is the integral of the largest principal strain rate. Thomason [74] used the above 

expressions in the derivation of the void coalescence condition. 

For an array of void nucleating particles of diameter Dp and spacing dg, setting the initial 

void radius as I?p/2, eqns. (1.12) and (1.13) leads to the following expression for fracture 

strain: 

In {dgJDp) (1 - n) 


£/ = 


sinh [(1 - n)((Ta + <7b)/(2cr/v'3)] 
In (dg/Dp) 


[McClintock, 1968] 
[RiceandTracey, 1969] 


(1.18) 


(1.19) 


0.28 exp(1.5crjfejfc/a) 

Some research workers have used the above expressions for the prediction of ductile fracture 
with the help of experimentally determined parameters. The major limitations of this ap- 
proach are that it ignores the effects of void nucleation (as the expressions are derived from 
a pre-existing finite size void) and void coalescence. Thus, e/ is usually overestimated. 


Void Coalescence Model 

Thomason [74] modeled void coalescence phenomenon as plastic instability due to necking of 
the inter void matrix. The sufficient condition for plastic instability of the inter-void matrix, 
as given by him, is 

Cn-An - (Ti = 0 (1.20) 

where An is the area fraction of the inter-void matrix perpendicular to the direction of maxi- 
mum principal stress ai and cr„ is the plastic constraint stress. He considered a geometrically 
equivalent square prismatic void and used the upper bound method to obtain an expression 
for the plastic constraint stress in terms of the void dimensions, inter-void spacing and the 
yield stress. 



Thomason’s Fracture Criterion 

Thomason [74] combined the results of Goods and Brown [39] on void nucleation, Rice and 
Tracey [41] on void growth and his own on void coalescence to arrive at a fracture criterion. 

He used expressions (1.14 to 1.16) for the void dimensions to express the void coalescence 
condition (eqn. 1.19) in terms of the void growth strain ef and the hydrostatic stress am- 
By superposing this condition on the void nucleation relation (eqn. 1.10), he obtained the 
fracture criterion as a graph of fracture strain (e/ = e” + ef) versus the triaxiality (i.e. the 
ratio of hydrostatic stress to the generalized stress). 

The main drawback of this approach is the use of expressions (1.14 to 1.16) for void 
growth. These expressions are based on integration procedure which assumes that the prin- 
cipal directions of strain rate remain fixed throughout the strain path. This is, in general, 
true only for the case of small strain and rotation. As a result, this approach cannot be used 
when the strains and/or rotations are large. 

Continuum Damage Mechanics Models 

In these models, material behavior is represented by a plastic potential which includes damage 
as an internal variable. The damage variable quantifies the intensity of voids, which can be 
identified as either the void volume fraction or the area void fraction. The basis of continuum 
damage mechanics models rests on the theory of continuum thermodynamics. In this section 
some of the continuum damage mechanics models are discussed. 

Lemaitre [47] proposed a continuum damage mechanics model for void growth in elasto- 
plastic materials with area void fraction as the damage variable using the concepts of effective 
stress and strain equivalence. His model is based on a simple damage growth law in which 
the damage rate is linearly dependent on the plastic strain and the thermodynamic force 
corresponding to the damage. It does not account for void nucleation. Note that, as far as 
void coalescence is concerned, it has to be incorporated as an additional condition in terms 
of the continuum parameters. This condition, which serves as a fracture criterion, has to : 

be based on an appropriate micro model. Thomason’s limit load model seems to be a good I 

candidate for this purpose. j 

More recently Dhar [46] and Dhar et al. [75] extended Lemaitre’s [47] continuum damage i 

t 

mechanics model to incorporate a damage growth law which accounts for both the void | 

nucleation as well as void growth. The law is based on the experimental results of Le roy et | 

al. [48]. They also modified Thomason’s [74] void coalescence condition to make it suitable 
for the case of large strain and/or rotation. They combined the extended continuum damage 
mechanics model and modified void coalescence condition to arrive at a criterion for prediction 
3 f micro crack initiation. The salient features of the above approach are discussed below. 

In this model, damage D is identified as an area void fraction at a point in a plane. Thus, 



it is defined bv 


D = 


AA, 

I^A 


( 1 . 21 ) 


where Aj 4 is the infinitesimal area around the point in some plane and AAt, is the area of 
the voids in the plane containing £\A. The conjugate variable corresponding to D is -Y, the 
rate at which the elastic energy releases during the damage growth at constant stress. For an 
isotropic material, -Y is given by [47] 


-y = 


2E{l-Df 



( 1 . 22 ) 


where 

/(^) = + 3(1 _ 2u) (1.23) 

a 6 a 

Here, E is the Young's modulus, u is the Poisson’s ratio and (tTm/^) is called the triaxiality. 
The damage growth law is given by 


b = CdI + {ai + asD) (-Y) t 


(1.24) 


where oi, as and Cd are material constants. Here the first term, which is independent of 
— y, represents the damage evolution due to void nucleation while the second term, which 
is dependent on -y, represent the evolution of damage due to void growth. For SAE 1090 
steel, Dhar et al. [75] determined the constants fitting eqn. (1.23) through the experimental 
results of Le Roy et al. [48]. They obtained the following values: 

Cd = 1.898 X 10-2 

Oi = 9.8 X 10-^ MPa-^ 
as = 1.84 MPa-i 

While getting these constants, they used Bridgman’s relation [49] to express the triaxiality as 
a function of strain. 

Dhar et al. [75] used finite strain expressions for the void dimensions and inter- void spacing 
to express Thomason’s void coalescence criterion (eqn. 1.19) in the following form: 




0.1 -f 


1.2 


{1 — exp (— e/2)}*^-® 


exp (-e) cr = 0 


(1.25) 


where e, the generalized strain, is the integral of i along the path lines. To find the crit- 
ical value of damage parameter (Dc), they applied this criterion to various geometries and 
loading conditions using the elasto-plastic finite element analysis. They observed that Dc is 
independent of geometry or loading and hence that can be used as a material property for 
the prediction of micro-crack initiation. They reported that the value of Dc is 0.05 for SAE 
1090 steel. 



Empirical and Semi-Empirical Models 

In the absence of reliable quantitative models for incorporating the phenomena of void nu- 
cleation, growth and coalescence in materials undergoing large plastic deformation in metal 
forming processes, many empirical relations based on some phenomenological models have 
been proposed. This section reviews some of them. 

Freudenthal [50] postulated that generalized plastic work per unit volume is the critical 
parameter and is expressed as 

f adl = C\ (1.26) 

Jo 

Cockcroft and Latham [51] modified the generalized plastic work criterion to take care of 
the effects of change in the neck geometry which is observed in tensile tests. The criterion is 
given by 

J — j aide — C 2 (1-27) 

Here ai is the maximum normal stress. 

Oh et al. [19] modified Cockcroft and Latham criterion [51] by considering the ratio of the 
maximum tensile stress to the generalized stress. This criterion is expressed as 

r ^de = Cs (1.28) 

Jo O' 

The main drawback of the above three criteria is that the effect of hydrostatic stress is 
not considered. Osakada and Mori [52] developed a ductile fracture criterion which includes 
the effects of both plastic strain and hydrostatic stress and expressed it as 


< o + 6 + hafn ^ de 


(1.29) 


< • >= • 


< • >= 0 • < 0 


where a and b are material constants. 

Norris et al. [ 53 ] developed a fracture criterion based on their experimental results and 
finite difference results of several test geometries and expressed it as 


/o (1 - Cam) 


(1.30) 


where C is a constant. 

In the above equations Ci, • • • C 5 are the critical values accumulated over the strain path 
to fracture. 



1.2.3 Defects in Forging 

Different type of defects are observed in industrial metal forming processes which are tabulated 
by Johnson and his co-workers [59, 64, 66]. The limit of forming processes is governed by 
the material properties, friction, temperature, rate of deformation and geometry. It is a 
well established fact that the tensile triaxiality is the main factor contributing to the crack 
initiation in those processes. 

There are several defects found in open die forging namely barreling, tears or surface 
cracks, centre burst, thermal cracks, end grains etc. Barreling occurs due to friction between 
the part and the die or when there is a large temperature difference between the part and the 
die. This leads to greater deformation in the midsection of the part than in the constrained 
ends [58]. The amount of barreling is more for unlubricated conditions [56]. Thermal cracks 
are the result of non-uniform temperatures within a metal. When flow lines reach a surface 
perpendicularly, exposing the grain boundaries directly to the environment, this condition is 
known as end grains. 

Cracks due to tangential velocity discontinuities or edge cracks (i.e. at the work-die 
interface) may arise in the material which is insufficiently ductile or when friction at the 
interface is high. In simple upsetting, longitudinal cracking occurs at the barreled surface due 
to presence of secondary tensions. A central cavity is a rupture occurring at the centre of the 
billet [54]. 

Thomason’s [55] uniaxial compression test results support the hypothesis that a metal 
must reach a state of tensile plastic instability before ductile fracture. He observed ductile 
fracture occurred on the surface of uniaxial compression specimens without any evidence of 
external local necking. 

Kobayashi [56] conducted a number of experiments with different height-diameter ratios 
in varying friction conditions. He reported that normal fracture occurred when friction was 
severe and the axial stress was tensile and when friction was less severe, the axial stress at 
fracture was compression and the fracture was of a shear type. 

Kuhn and Lee [57] conducted experiments on steel cylinders under varying conditions. 

A refined measurements of strains revealed that the axial strain suddenly reached a plateau I 

and remained constraint while the hoop strain continued to increase. The axial strain then i 

[ 

resumed its increase and did so at the same rate as the hoop strain. They interpreted this | 

strain perturbation as a localized instability in strain and hence localized necking of material, | 

[ 

between inhomogeneities, which eventually led to fracture. Also they metallographically | 
detected void formations prior to fracture. 

Sowerby, Dung et al. [60,62] examined the capability of McClintock’s void growth model 
, fp predict damage accumulation in the upsetting of steel specimens using rigid-plastic formu- 
llation and found it appropriate for homogeneous metals. Predeleanu et al. [63] conducted 



experiments and computer simulations using their damage criterion and found that material 
damage reaches a critical value inside the specimen and not at the free surface. They also 
detected voids and micro-cracks within the interior of the specimen before surface fracture. 

Clift et al. [61,65] used 3-D elasto-plastic finite element analysis using 8-noded brick 
element. They modeled boundary friction conditions by defining an extra layer of elements, 
referred to as the friction layer on the required surfaces. But they didn’t report effect of friction 
in their study. The numerical results were used to give predictions of fracture initiation site for 
a total of nine fracture criterion. They compared these predicted sites with experiments and 
showed that generalized plastic work per unit volume criterion is capable of estimating the 
experimental fracture initiation sites for the complete range of processes considered. They 
re-confirmed that the fracture criterion proposed by Ghosh was not successful in any one 
of the cases examined. Later Zhu et al. [67,69] also disapproved Ghosh criterion in their 
simulation on side pressing of cylindrical billets. The main reason for the failure of Ghosh’s 
criterion is that it was proposed for sheet metal forming, not for bulk metal forming. 

Lin and Lin [68] used strain energy density failure criterion in conjunction with coupled 
thermo-elasto-plastic large deformation model to analyze ductile fracture in upsetting of the 
cylindrical specimen under various geometrical and frictional conditions. They reported in- 
creased damage with higher friction. 

Semiatin et al. [72] conducted hot forging trials on Ti-6A1-4V specimens and validated the 
results with finite element analysis and suggested that the Rice and Tracey void growth model 
which highlights that the influence of hydrostatic stress and Cockroft and Latham criterion 
provided good estimate of surface fracture. 

Atkins [70] found that well lubricated cylinders developed inclined cracks found in shear 
while unlubricated ones showed vertical cracks formed in tension. He also mentioned that the ! 
greatest strains to fracture occur at the least values of tensile hydrostatic stress. > 

Recently, Kim, Im and Geiger [73] investigated the ductile fracture criterion based on I 
work hypothesis and Cockcroft and Latham criterion, and carried out experiments and rigid- | 
visco-plastic finite element analysis studies on simple upsetting of aluminium billets. They 
deduced that Cockcroft and Latham’s criterion gave a more reasonable prediction for crack 
initiation site than work hypothesis and the former was useful in processes like cold forging 
in which the influence of the induced circumferential tensile stress on failure is dominant. 

Gouveia et al. [71] used finite element analysis and experimental results to analyze four 
ductile fracture criteria in forging. They reported that Preundenthal criterion predicted much 
higher values of damage inside the specimen and thus located fracture initiation site at the 
center of the billet. The drawback of Cockroft-Latham criterion, they suggested, was if the 
largest principal stress aj is temporarily negative led to accumulation of negative damage. 
They concluded that the Oyane’s criterion which is formulated on the basis of the void growth 
model and the theory of the plasticity of the porous media, was the best among the criteria ; 



considered. 

The above studies show that published ductile fracture criteria have had only limited 
success in predicting the fracture in forging process. Further, different experimental set-ups 
seem to give different values of the critical material parameters required in the criterion. A 
universal criterion which can be applied to all forging problems has yet to come up. 

1.3 Scope and Objective of the Present Work 

In the present work, damage in axisymmetric cold forging process is studied. An axisymmetric 
large deformation elasto-plastic finite element code is developed. The Updated Lagrangian 
formulation, which is convenient for handling material and geometric non-linearities, is used. 
New incremental objective stress measure and logarithmic strain measure are used instead 
of Jaumann stress rate and Green-Lagrange nonlinear strain tensors used in most of the 
literature. Forging problem is a displacement-controlled problem. Therefore, to accelerate 
the convergence of iterative scheme, arc length method is used in conjunction with modified 
Newton-Raphson iteration technique. The material is assumed to be elasto-plastic strain 
hardening and yielding according to von Mises criterion. Inertial forces and body forces are 
not included in the present work. Coulomb’s law is used to model friction at the die-workpiece 
interface. 

Review of the published research reveals that there have been several attempts to predict 
the occurrence of cracks in forging. But most of the fracture criteria used are empirical ones 
and they have had only a limited success. Ductile fracture in metal forming processes is 
now established as a void nucleation, growth and coalescence phenomenon. Dhar et al. [75] 
combined the continuum damage mechanics model of Lemaitre [47] with the void coalescence 
criterion of Thomason [74] to propose a criterion for micro-crack initiation in terms of the 
critical value of the damage parameter. This criterion incorporates all the features of mi- 
croscopic description of ductile fracture, namely void nucleation, growth and coalescence. In 
the present study, this criterion is taken to predict damage distribution in AISI-1090 steel 
workpiece. This material is chosen because of availability of the values of constants used in 
the criterion. 

The code is checked for mesh and step size convergence. Next, a detailed parametric 
study of damage distribution is carried out to show the effect of coefficient of friction at 
the die-workpiece interface and the height to diameter ratio of blank. For typical values of 
coefficient of friction and the height to diameter ratio, axial, circumferential and hydrostatic 
stress distributions are studied. 



1.4 Organization of the Thesis 


The thesis is organized as follows. In the second chapter, the mathematical modelling of 
forging process, the axisymmetric finite element formulation and the boundary conditions are 
presented. The damage formulation is then discussed. Numerical scheme is presented at the 
end of chapter. In the third chapter, results for a typical case are presented. It also includes 
discussion on the parametric studies. 



Chapter 2 


Mathematical Modeling and Finite 
Element Formulation of Axisymmetric 
Forging Process 


In this chapter, mathematical model and finite element formulation of a cold forging pro- 
cess are developed. The process is modelled as an axisymmetric problem. The description 
of motion, stress measures, strain measures and constitutive laws used in the formulation of 
the governing equations for static large deformation elasto-plastic problems are presented in 
sections 2. 1-2.7. The finite element formulation is discussed in section 2.8. The boundary 
conditions for axisymmetric forging process are discussed in section 2.9. The damage for- 
mulation is outlined in section 2.10 and the solution procedure for solving the finite element 
equations is presented in section 2.11. 

2.1 Updated Lagrangian Formulation 

In the study of the deformation of a body subjected to external loading, often the original 
undeformed and unstressed state of the body is used for the formulation of its equation of 
motion. This is known as Lagrangian formulation. This formulation is convenient for small 
deformation, which is the case in majority of engineering problems. In such cases, the de- 
formed configuration does not deviate much from the original one and hence the deformation 
can be described by an infinitesimal strain tensor, also known as the engineering strain tensor, 
for which the strain-displacement relations are linear. On the other hand, for large deforma- 
tion problems, one has to use a finite strain measure, which is expressed by a nonlinear 
strain-displacement relation. Furthermore, the equations of motion when expressed in the 
reference configuration depend on the deformation. Hence, for large deformation problems, 
the Lagrangian formulation proves to be cumbersome with the governing equations being 
difficult to solve. In such cases, one solves the problem using an incremental method known 
as Updated Lagrangian formulation. In this formulation, it is assumed that the states of 



stress and deformation of the body are known till the current configuration, say at time t. 
The main objective is then to determine the incremental deformation and stresses during an 
infinitesimal time step, At, i.e from time t to t+At. Here, the current configuration is used 
as the reference configuration for obtaining the incremental values. Unlike in the Lagrangian 
formulation, an incremental strain tensor is used. This methodology is particularly useful 
for elasto-plastic materials because the stress-strain relationship in such materials is usually 
expressed in an incremental fashion. 


2.2 Kinematics of Finite Deformation 


The relative deformation gradient tensor at time t + At is defined by 


jp 




J+At 


d% 


t J 


( 2 . 1 ) 


where and denote the position vectors of the particle at times t and t + At respec- 
tively. When is not singular, the polar decomposition theorem [76] allows a decompo- 

sition of the form 


=‘*‘■1 , ( 2 . 2 ) 

where is an orthogonal tensor representing the material rotation and is a 

positive definite symmetric tensor called right stretch tensor. The right stretch tensor can be 
diagonalized by the following transformation to obtain the principal stretches : 

=‘+i Oii , (2.3) 

where is orthogonal and 

( 

‘+"A, =<+" £11 . (No Sum) (2.4) | 

I 

The tensor and the transformation matrix are used in the development of a | 

stress measure while the incremental logarithmic strain definition is obtained from the 

2.3 Stress Measures 

It is essential that an objective stress measure be used in the incremental theory of constitutive 
modeling to account for the rigid body rotation that may accompany deformation. The 
Cauchy stress tensor, which has great physical significance, is not objective and hence cannot 
be used directly in a constitutive equation. There are numerous objective stress measures, 
each with particular advantages and disadvantages. Some of them are discussed below. 



One of the most commonly used objective stress tensor is the second Piola-Kirchoff stress 
tensor which can be related to the Cauchy stress tensor using the concept of 

equivalent work between two configurations [78]. 


f^Oij — 
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(2.5) 


where *'p and are the densities at time t and t+At respectively and denotes 

the derivative of *Xi with respect to The second Piola-Kirchoff stress tensor is energy 

conjugate with the Green-Lagrange strain tensor. This energy conjugate pair is used in the 
equation for predicting displacements in a predictor corrector solution procedure (see section 
2.7). 

Another commonly used objective stress rate measure is the Jaumann stress rate <t. It is 
related to the Cauchy rate & by 


' a% At =Vy At Oik CW.kAt) ajk CWikAt) , (2.6) 

where 

^WijAt = -t Auj^i) , ' (2.7) 

represents the component of the incremental spin tensor. Here, tAuij denotes the deriva- 
tive of the incremental displacement vector tAu with respect to ‘a;. 

As pointed out by Metzger and Dubey [79], it is important that the stress measure used be 
compatible with the constitutive equation in addition to being objective. An incremental ob- 
jective stress measure (as against objective stress rate measures) to be used in the generalized 
Hooke’s law and in the elasto-plastic constitutive equation is developed below. 

Two Cartesian reference frames (see Fig. 2.1) are used. 

1. The fixed frame. 

2. A material frame which rotates and translates along with the material particle. 

The material frame is defined so as to coincide with the principal axes of the right stretch 
tensor at time t so that the initially orthogonal axes do not get skewed at time t + At. 

The first step in the derivation is the transformation of the components of the Cauchy 
stress tensor at time t from the fixed frame to the material frame: 


(2.8) 

The increment in the components of the Cauchy stress tensor with respect to the material 
axes, is added to to obtain the stress components at time t + At with respect to 

the material frame of reference: 
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The fixed frame is denoted by subscript F, 
the material frame is denoted by subscript M. 

Figure 2.1: Fixed and material reference frames 

The final transformation from the frame at time t+Atto the fixed frame give the components 
of the Cauchy stress tensor at time t + At : 

=‘+^‘ Qii (2.10) 

The equations (2.8), (2.9) and (2.10) can be combined together to give a relation between 
the components of the Cauchy stress tensor at times t and t + At and the increment in the 
components with respect to the material frame: 

=‘*^l Q„ ‘+“0™ +■ Ac.") ‘+“0,, (2.11) 

The proposed stress and the logarithmic strain measures satisfy the requirement of objectivity 
and lead to a physically consistent application of the usual constitutive equation. 

2.4 Strain Measures 

The incremental small strain is given by 


tAsij — +t Auj^i) 


( 2 . 12 ) 


The incremental Green Lagrange strain is a non-linear function of the displacement 


tA&ij — -\~t Auj^i -l-t Auk^ifAuk^) 


(2.13) 


These two strain measures occur in the virtual work expression at time t + At and its 
transformation to time t respectively. The components of the Green Lagrange strain tensor 
are invariant under rigid body rotation of the material unlike the small strain components. 

The following are the disadvantages of using one of the above measures in a constitutive 
law: 

1. The solution obtained is dependent upon the size of the increment in the Updated 
Lagrangian formulation unless the increment size is sufficiently small. 

2. The components of the strain tensors do not tend to infinite values when the principal 
stretches tend to zero. Therefore, a constitutive law, which ensures that the appropriate 
Cauchy stress components tend to negative infinity (as is physically realistic), even though 
the strain components remain finite, should be used. This difficulty can be avoided by using a 
strain measure whose components become minus infinity when the principal stretches become 
zero. 

The logarithmic strain measure introduced by Dienes [77] is free of the above disadvan- 
tages. The principal incremental logarithmic strain components, which will be used in this 
work, are defined by 

.a 4 = in(‘+“Ai) ii, ■ (2.14) 

where the A, are the principal stretches defined by equation (2.4). The logarithmic strain 
has the following additional advantage in elasto-plastic analysis. A loading test involving 
elasto-plastic deformation followed by elastic unloading reveals that the slope of the elastic 
unloading line is the same as that of the initial elastic line only when the true stress and the 
logarithmic strain measures are used in a constitutive law [80]. 

2.5 Elastic Constitutive Equation 

A constitutive law must satisfy the principle of material frame indifference [76]. The incre- 
mental elastic constitutive equation used is the generalized Hooke’s law relating the increment 
in stress components with respect to the material frame and the elastic part ([As®^]) of the 
principal incremental logarithmic strain components: 

(2.16) 

The tensor [C^\ for the isotropic case is given by 

= + (2.16) 

where A and /i are Lame’s constants. 

In case anisotropic behavior is modeled, the tensor C® has to be evaluated with reference 
to the directions of the principal stretches. 



The stress and strain measures used in a constitutive equation need not necessarily be 
energy conjugate with each other. However if it is so, the predicted response in a predictor- 
corrector scheme will be closer to the actual response. 

2.6 Elasto-Plastic Constitutive Equation 

As stresses developed in a material exceed the yield stress, the elastic constitutive relationship 
between the stress and strain tensors no longer remains valid. We now develop a relationship 
between the stress and strain based on the von Mises yield criterion and isotropic hardening. 
For an isotropically hardening material, the plastic potential is given by [81] 


F{(Jij,p) = aeq{c^ij) - CTyip) (2.17) 

note that 

F = 0 (2.18) 

represents the yield criterion. The plastic potential F depends on the Cauchy stress tensor 
aij through it’s second invariant aeg called as equivalent stress and defined by 

(2-19) 

where a ij is the deviate ric part of cry. Further, F depends on the variable yield stress of the 
material, ffy, through a hardening variable p. For the case of strain hardening, p is identified 
as the equivalent plastic strain and hence defined as: ^ 

P = £lg = J <Hq (2-20) I 

where ^ 

< = (|<i4 de?,) ’ ' (2.21) ; 

I 

Here, dsy is the plastic part of the incremental linear strain tensor dsy and the integration | 
in equation (2.20) is to be carried along the particle path. The dependence of ay on p (or | 
£^g) is normally approximated by a power law type of relationship 

a, - = K(4,r (2.22) 

Here, {ay)o is the yield stress at zero plastic strain, K is called the hardening coeflicient and 
n is called as the hardening exponent. 

The plastic part of incremental linear strain tensor (ds?- ) is obtained from the plastic 
potential using the following relation: 



where dX is a scalar. This equation is called the flow rule. Differentiation of equation (2.17) 
with respect to (7^ gives 


(9<7„ 2creg‘^*^ 

using this, one can determine dX as: 


(2.24) 


dX = dsl^ 


(2.25) 


Further, the hardening relationship and the yield condition can be used to express dX as: 


where 
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(2.26) 


= ^ = Kn(el,T-' (2.27) 

is the slope of the hardening curve. Substitution of equation (2.24) and (2.26) in equation 
(2.23) leads to the following constitutive equation 


p _ ^ ' 


(2.28) 


This constitutive relationship between the deviatoric stress tensor and the plastic part of linear 
strain tensor is not really convenient for the Updated Lagrangian formulation for which the 
incremental stress-strain relationship is needed. This can be obtained from equation (2.28) 
as follows: 


, p 3 d<ye.q j 

Note that, from equations (2.17) and (2.24), we get 
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Substitution of equation (2.30) in equation (2.29) leads to the following incremental plastic 
stress strain relationship: 

9nr'. ./t! . 

(2.31) 
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The incremental elastic stress strain relationship (equations (2.15)-(2.16)) can now be written 



— '^[—^dakk^ij + (1 + u)daij] (2.32) 

where def^ is the elastic part of d£ij, E is the Young’s modulus and u is the Poisson’s ratio. 
Adding the two relationships, we get 


d£^j = de^j + de^i^ 


(2.33) 
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This is the incremental elasto-plastic stress strain relationship needed in the Updated La- 
grangian formulation. However, it is the following inverse relationship which is more useful: 


where 


dffij = dski 


(2.35) 
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and IX is one of the Lame’s constant (also called shear modulus) appearing in equation (2.16). 

Note that the stress increment appearing in equation (2.35) must be an objective stress 
increment in the sense that daij must reduce to a zero tensor in the event of the increment 
being a pure rigid rotation. The incremental objective stress measure to be used in the present 
work has been described in section 2.3. 

The relationship (2.35) has been derived assuming the increment size to be small and 
suing the incremental linear strain tensor as the strain measure. When a large size increment 
is to be used along with the incremental logarithmic strain measure (defined in section 2.4), i 
the derivation is similar. The relationship (2.35), when the stress and strain measures are ! 

replaced respectively by the incremental objective stress measure (of section 2.3) and the ; 

incremental logarithmic strain measure (of section 2.4), takes the following form with respect \ 
to the material frame. ! 


M 
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(2.37) 


where 


2(3/1 + irn(>4,)”-i)V?, 

Here, H has been replaced by the expression (2.27) and the left superscript t has been added 
to make it explicit that these quantities are to be evaluated at that time. 




(2.38) 



2.7 Incremental Updated Lagrangian Formulation Carte- 
sian 

The objective of the Updated Lagrangian formulation is to establish static equilibrium in the 
configuration at time t + At when all static variables at time t are known.The principle of 
virtual work requires that 
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where 

= Cartesian component of the Cauchy stress tensor at time t + At, 


rt+At _ 1 ( 

0 £ij — 2 af+Atj.^ T gt+Aix^ J ; 


= Cartesian component of the virtual displacement at time t + At, 

= Cartesian coordinate of a material point at time t + At, 

e+Aiy = volume at- time t + At and 

r ^2.40) 

Jt+^tSf 

where 

t+£^tf. = Cartesian component of external traction at time t + At 
= Surface at time t + At with traction specified. 

The main difficulty in the application of equation (2.39) is that the configuration at time 
t + At is unknown. An elegant way of formulating the problem is given by Bathe [82]. The 
virtual work expression at time t + At is transformed to an integral over the volume at time 
t by using the second Piola-Kirchoff stress and Green Lagrange strain measures, which are 
energy conjugate to each other, and the principle of conservation of mass, it is assumed that 
the external load term (2.40) is deformation independent for the formulation of the governing 
equation. The expression (2.39) after the transformation becomes 
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The Green-Lagrange strain tensor is defined by 
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where denotes the derivative of with respect to ^Xj. 

The second Piola-Kirchoff stress tensor can be decomposed as 


=; Sy +, AS,, =• < 7 ,, +, AS,, 


(2.43) 



Since 


= S{\ +( Awi) = S{tAui) 

the virtual Green-Lagrange strain tensor can be decomposed as 


— 5{iAeij) — 5{tAsij +j ^Vtj) 


where 


t^Vij ~ “ht 


Therefore, equation (2.41) can be written with incremental decomposition as 
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(2.45) 

(2.46) 
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(2.48) 


The above equation is linearized by neglecting the second integral, which is a higher order 
term, and approximating jAS'y as As^P 
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The linearized equation, when solved, will yield only approximate displacement, strain 
and stress fields. The approximate quantities are denoted by a right superscript (1). The 
error to the approximation involved is calculated from equation (2.39) as 

Error f (2-50) 

Jt+AtvW 

This error is generally minimized by an iterative predictor-corrector scheme described in 
section 2.11.2. 

The next section discusses the finite element formulation of the governing equation (2.49) 



2.8 Finite Element Formulation 

2.8.1 Matrix Notation 

Matrix notation is used in the development of the finite element equations for the convenience 
of computer implementation. 

For the axisymmetric case, the components of the tensors, jAe, jAt? and tAe^ are repre- 
sented in array form as follows: 



(2.51) 

t{^V} — )t AliriZ) t^'^z,ri t^'^z,zi } 

(2.52) 

t{AeY'^ = {tA4,tA£,\jA4_eA4j 

(2.53) 


The components of all stress tensors are written as arrays with respective components placed 
in the following manner: 


[aY = {cTrn cTrz, (^eo} (2.54) 

Two matrix forms of the tensor (given by equation (2.38)) are needed: 

1. for the constitutive relationship : t{ A5} = {{As}, 

2. for the constitutive relationship : t{ Acr^} = ‘[(7^^] t{A£^}. 

They are given by 





tH +t {aY[C^Y{a} J 

In the above equations. 


(2.55) 

(2.56) 


‘{4 = jjAV'} (2.57) 

and the array ^{a} represents the deviatoric part of the Cauchy stress at time t. 

For an isotropic material, the matrices [C^ ] and [C^] for the axisymmetric case are given 
by 
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Equation (2.49) can be written in the following form owing to the symmetries in 
ffAe}, i{A77} and : 
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where the matrix *[T] is given by 
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2.8.2 Finite Element Equations 


(2.60) 


(2.61) 


The domain is discretii^ed into a number of elements and the incremental displacement field 
is approximated over each element by 


,{Au} = I } =■ [4],{A,<}' (2.62) 

where the element incremental displacement vector t{A'u}® for an n-noded element is given 

by 


= {tAul , t^ul , A< A<} (2.63) 

The quantities ^A^^^^ tAu% stand for the unknown incremental displacements of node i 
in r and 2 ; directions respectively and the matrix *[$] is defined by 


where 
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= 0, ‘«2, 0 w„o} 


‘{$,}’' = {0,‘iVx,0, % 0 , ‘JV„} 
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The *Nt, which are functions of {^r,^z) are called shape functions. The number of nodes 
and the corresponding shape functions are chosen on the basis of convergence criteria. For 
the problem under consideration, 8-noded serendipity element is used. 

The strain field is expressed in terms of the nodal displacements by differentiating (2.62) 
and using the expression (2.46) and (2.47). This leads to 

t{Ae} [BlUAuY (2.66) 

t{Ari} =‘ (2.67) 

where 

■ ■ 

W = ■ ( 2 - 68 ) 

and 

■ ■ 

= ’ (2.69) 

Using equations (2.62), (2.66) and (2.67), the contribution to the integral (2.60) over a 
typical element e with volume U® is 

{ (.{Au}*’’) U ‘[BLT'[C^'’']‘[BL\<fV''^ t{A«}'+ 

5 (.{Ati}'’’) U ‘VBxf ,{A«}'+ 

S (,{A«}'’') ‘H cfV’'^ = S (.{Au}'’') ‘+“{f'}‘ (2.70) 

The contribution to the term is expressed in terms of the elemental external force 

vector using a standard procedure. Since the variation in the displacement vector 

is arbitrary, the above equation can be written as 
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where the elemental stiffness matrix ‘[F]® is given by 


‘IK]‘ =■ IK,Y +‘ [K„i\‘ 

(2.72) 
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(2.73) 

\KmlY= [ 

JiV‘ 

(2.74) 

and the elemental internal force vector is 


'{fY= [ 

Jtve 

(2.75) 

The elemental stiffness matrix ‘[F]® along with the elemental force vector *{fY 
i+AtijT’j.e assembled to obtain the global equation. 

‘[JfHAu} +‘ {/} =‘+^‘ {F} 

(2.76) 

Decomposing equation (2.76) can be written as 


^[K]t{Au} +‘ {/} =‘ {F} +t {AF} 

(2.77) 

Here, ‘{F} is the (global) external force vector at time t and t{AF} is the (global) incre- 
mental force vector from time ttot + At. In Updated Lagrangian formulation, it is assumed 
that the equilibrium equations are satisfied exactly at time t. Thus 

■{/} =‘ {F} 

(2.78) 

‘[A-],{A«} =, {AF} 

(2.79) 


This equation is solved to obtain the incremental displacement vector t {Au}. 


2.8.3 Static condensation scheme 

When the Gauss elimination method is applied to the problems consisting of contact boundary 
conditions after incorporating them, there is a possibility that the pivot may become zero. 
Storage of the coefficient matrix in the full form is also not desirable. Since the contact 
condition needs to be applied in an iterative manner, solving the full storage equations in 
all the iterations takes an enormously large amount of time, especially when the number 
of degrees of freedom is high. To reduce the solution time, it is desirable to condense the 



coefficient matrix and the right side vector to the size involving only those degree of freedom 
to which the contact boundary condition is to be applied. This condensation procedure is 
described in the next paragraph. 

First the total nodes arc divided into two categories. The nodes at which contact condition 
is to be applied are called the nodes of type 2. The rest of the nodes are called as nodes of type 
1. To facilitate the condensation, the nodes are numbered in such a manner that the nodes 
of type 1 are numbered first and then the nodes of type 2 are numbered. Then the essential 
boundary conditions (i.e, the boundary conditions on the incremental displacement vector, 
which is the primary variable in this case) for the nodes of type 1 are applied to equation 
(2.79). The essential boundary conditions of the forging process are described in section 2.9. 

Next, these equations are partitioned as follows^: 


[Ku] [K.,] ] f {AuJ 1 _ f {AF.} 1 
[if2ij [KJi J \ {Atij} / - \ {Af,} / 
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where {Aui} denotes the incremental displacement vector corresponding to the nodes of 
type 1 and {AU 2 } contains the incremental displacements of the nodes of type 2. Equation 
(2.80) can be separated as follows; 


[ii-,.]{Aui} + [/s:i2l{A«2} = {AF,}, 

(2.81) 

[ifjiHAu,} + lKi2]{Ma} = {AFj}. 

(2.82) 

solving for {Au,} from equation (2.81), we get 


{Aui} = [jru]“^({AFi} — [ii'i2]{Au2}). 

(2.83) 

Substitution of this expression for {Aui} in equation (2.82) and rearrangement of the 
resulting equation leads to the condensed set of approximate equilibrium equations: 

m{Au} = {AF} 

(2.84) 

where 
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{AF} = {AF 2 } - (ifjiKBl, 

(2.87) 1 


^In this section, the superscripts/subscripts denoting time have been omitted for the sake of convenience. 
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{B} = (ii-„l-‘{An}, (2.89) 

The number of equations in the condensed set (2.84) is equal to the number of degrees 
of freedom of type 2. The coefficient matrix [ilT] and the right side vector {jP} are evaluated 
from equations (2.85-2.89) using the partitioned matrices and vectors of equation (2.80). Note 
that while evaluating [/I] and [B], it is not necessary to invert the matrix [Ku]. Instead, one 
can solve the following equations by the Gauss elimination method: 




(2.90) 


[Ku]{B} = {ABi}, (2.91) 

To reduce the computational time further, one can store the upper triangular from of [ATn]. 
In this way, one can obtain columns of [A] by performing the Gauss elimination and back 
substitution operations on the corresponding columns of [^^ 12 ]- The vector {B} is obtained 
in similar fashion by performing the Gauss elimination and back substitution operations on 
the vector {ABj}. 

2.8.4 Determination of Stresses 

The evaluation of the stress components (at the Gauss points of the elements) is done by the 
following stepwise procedure: 


1. Calculate the relative deformation gradient 


‘+“[F] = 
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It should be noted that the position vector ‘cc corresponds to the equilibrium position. 


2. Decompose the relative deformation gradient as and determine 

and using equation (2.3). 

3. Determine the principal incremental logarithmic strain t{As^} using equation (2.14). 



4. Calculate e{A(7^^}using equation (2.37). The integration of the constitutive equation 
is performed using the Euler forward technique described below. 

5. Use equation (2.11) to calculate the Cauchy stress components at time t + At. 

2.8.5 Integration of the Constitutive Equation 

Different techniques exist for the integration of the constitutive equation. A simple and robust 
technique is the Euler forward integration scheme described below. 

Assume that the principal incremental strain components have been calculated and the 
state of the Gauss point at time t (elastic or plastic) is known. 

•If the state at time t is elastic, 


1. Calculate the stress increment assuming elastic behavior 

= [C®] ,{A£q. 

2. Calculate using (2.11). 

3. Determine ^aeg and using (2.19). 

4 elastic behavior holds, or if =^cy^q the Gauss point is 

neutrally loaded. Return. 

5. If I a transition from elastic to plastic has occurred. Calculate the 


Ratio = 


{*+^‘cre, -‘ffe?) ’ 


Change the state to plastic and use the sub-incrementation method. In this method, the 
stress components with respect to the material frame are updated after each sub increment by 
the increment in stress components corresponding to the elasto-plastic strain sub-increment. 

The [C^^] corresponding to the last updated state is used: 


where 



Radius 


upper; PLATEN 
A: 


% 

<u 

X: 


D;.. 


lower; PLATEN 


Figure 2.2: Domain of the Problem 


HAi{^M}(o) ^ + Ratio[C^]t {Ae^}, 

i+At{^£:p}(o) evaluated at 

Use equation (2.10) to find Return. 

• If the state at time t is plastic, the sub increment method described in (5) above is 
applied with Ratio set to zero. 

This describes the Euler forward scheme for the elasto-plastic model described in section 

2.6 

2.9 Boundary Conditions for Axisymmetric Forging of 
Cylindrical Disc 

Figure 2.2 shows the domain of the problem. Due to symmetry, only the upper half of the 
workpiece is used for analysis. The boundary conditions for the domain of this problem are 
discussed below. 



2.9.1 The Workpiece-Die Interface (AB) 

As stated in section 1.2 of chapter 1, the friction at the workpiece die interface is modelled 
by the Coulomb’s law. Therefore, a node of the workpiece at the interface either sticks to the 
die or slips relative to the die depending on the magnitude of the stress in r-direction and the 
coefficient of friction /i. Therefore, the boundary condition in r-direction can be stated as: 

- IJL - 0 i/ > fx (2.93) 

or 

tAur = 0 if < /X (2.94) 

Since both U and are negative, the first condition can be stated as 

= 0 if (2.95) 

Since nodal forces are proportional to the stress components, this condition can be ex- 
pressed in terms of the nodal forces also. In fact, this is the form, which is convenient in finite 
element formulation. Thus, the boundary condition in r-direction becomes: 

Natural boundary condition: 

= 0 if > fx (2.96) 

or 

Essential boundary condition: 

tAur = 0 if < /X (2.97) 

Since this is a displacement control problem, z-component of incremental displacement 
vector is specified at the interface. Thus, the boundary condition in z-direction is: 

Essential boundary condition: 


(Axxz = specified (2.98) 

2.9.2 The Free Surfaces (BC) 

These boundaries are traction free surfaces. Therefore, the boundary conditions on BC are. 
Natural boundary conditions: 

AFr = 0, xAFz = 0. 

Here, tAFr and tAF^ are the components of the incremental force vector. 


(2.99) 



2.9.3 The Plane of Symmetry (CD) 

Because of symmetry, 2-component of incremental displacement vector and r-component of 
incremental traction/force vector are zero. Thus the boundary conditions on CD are: 
Natural boundary condition: 


AFr = 0 . ( 2 . 100 ) 

Essential boundary condition: 

Au^ = 0 . ( 2 . 101 ) 

2.9.4 The Axis of Symmetry (DA) 

Again, due to symmetry, r-component of incremental displacement vector and z-component 
of incremental traction/force vector are zero. 

Natural boundary condition: 


AF, = 0 . ( 2 . 102 ) 

Essential boundary condition: 

Aur = 0. (2.103) 

2.10 Damage Formulation 

A continuum damage mechanics criterion proposed by Dhar [46,75] is used to predict micro- 
crack initiation in the present work. This has been discussed in section 1.2.2. Implementation 
of this criterion is discussed below. 

2.10.1 Implementation of Dhar ’s [46,75] Model 

In this model, damage D is identified as an area void fraction at a point in a plane. It is 
defined by equation (1.21). 

The damage growth law which gives rate of change of D is given by equation (1.24), which 
is reproduced below: 


t) = CA + (oi + CL 2 D) (-F) k (2.104) 

where -F is the rate at which the elastic energy releases during the damage growth at 
constant stress. For an isotropic material, it is given by [47] 



(2.105) 


where 


-V = 


2B(1-D)2 



/(^) = + 3(1 _ 2u) (2.106) 

ff 3 O' 

As already described in section 1.2.2, the constants fitting equation (2.104) are found for 
SAE 1090 steel as: 

Cd = 1-898 X 10“^, ai = 9.8 x 10~'‘ MPa~^ and oa = 1-84 MPa"^ 

Equation (2.104) can be written to give the value of increment in I? at a point in the 
domain as 


tAU = CdA£ + (ai +a2‘D)(-‘+'^*F)tA£ 

where 

t+Atpj^ t+At„ 

_(t+A«y^ ^ ^ ff 

^ ’ 2E [I - Wf ^ ' 

where function / is as defined in equation (2.106). 

Since, the strains and stresses are calculated at the Gauss points, value of LD can be 
found at every Gauss point in an element. This can then be used to determine damage 
growth variable D at time t + At. 


(2.107) 


(2.108) 


t+Atj^ = ‘D + jAD (2.109) 

The value of D is taken as zero at the beginning. Thus, a damage distribution in the 
domain can be obtained. The crack initiates wherever the value of D reaches the critical limit 
of 0.05 (for SAE 1090 steel). 

2.11 Solution Procedure 

Equation (2.84) represents only an approximate equilibrium equation at time t + At { the 
approxims^tion is mostly du6 to tho linsarisation and simplification involved in the steps 
between equations (2.48) and (2.49). A solution of such an approximate equation may involve 
a significant amount of error and depending on incremental displacement step, it may become 
unstable. Therefore, it is necessary to modify equation (2,84) to turn it into iterative problem 
capable of providing a solution with desirable accuracy. Amongst various iterative techniques 
[78], the modified Newton-Raphson algorithm is the most effective as it offers fast convergence 
with less computation. This algorithm is described in section (section 2.11.2). 



As stated in section 1.2 of chapter 1, in the present problem, the friction at the workpiece- 
die interface is modelled by the Coulomb friction law. The corresponding friction condition 
is given by equations (2.96-2.97), which can be applied only in an iterative fashion. Thus 
there are two sets of iterations in the solution procedure of the present problem. First, for the 
specified incremental force/ displacement, the friction iterations are carried out to determine 
the status (sticking or slipping) of the interface nodes. This algorithm is described in section 
2.11.1. Then, corresponding to this status, modified Newton-Raphson iterations are carried 
out to minimize the error arising out of linearisation and simplification of the equilibrium 
equation. 

Forging is a displacement control problem i.e the deformation is controlled by a specified 
displacement (at the workpiece-die interface) rather than by a specified traction. In such 
problems the rate of convergence of the modified Newton-Raphson iterative scheme is usually 
slow. To accelerate the rate of convergence, the arc length method is used in conjunction 
with the modified Newton-Raphson scheme. This scheme is explained in section 2.11.3. 

Finally, the numerical integration scheme for the elemental coefficient matrices and right 
side vectors and some divergence-handling procedures for the modified Newton-Raphson 
scheme are presented in sections 2.11.4 and 2.11.5 respectively. 

2.11.1 Friction Algorithm 

This section describes the iterative procedure to apply the friction boundary conditions (eq. 
(2.96-2.97)) to the linearized and condensed set of finite element equations (2.84). The itera- 
tive algorithm is as follows: 

(I) First step: 

1. In the first iteration, it is assumed that all the interface nodes are sticking. Thus, the 
boundary condition given by the equation (2.97) is applied to all nodes. 

2. Next these equations are solved to obtain {tAu 2 }- The vector {tAui} is determined 
from equation (2.83). Thus the whole incremental displacement vector {tAt6} is known. 

3. Then, the reactions are found by multiplying the original coefficient matrix with the 
incremental displacement vector. 

4. At the end of the contact iteration, the slipping nodes at the interface are identified 
using the condition 



(II) Second step: 

1. Now the condition (2.96) is applied to all the slipping nodes. This is done as follows: 

If i is the type 2 node number of a slipping node, then application of condition (2.96) to 

equation (2.84) means replacing {2i — 1)‘^ row of ‘[iiT] by the combination: 

^K 2 i-i,j - for j=l, number of degrees of freedom of type 2 

and setting (2i - row of ‘[AF] to zero. 

2. Then the resulting equations are solved to find {tA'U 2 }. The vector {fAui} is found as 
before. 

3. Now the reactions are found by the method described above. 

4. Finally, a check is made to find whether any additional nodes are slipping by using the 
condition 

(III) Further steps: 

The second step is repeated till the condition < fj, is satisfied at all the 

nodes. 

2.11.2 Modified Newton-Raphson Scheme 

The condensed set of approximate equilibrium equations (2.84), when the superscript/superscript 
denoting the time are restored, takes the form: 

‘[F]dA«} =t {AF} (2.110) 

The modified Newton-Raphson iterative scheme is used to convert this problem into an 
iterative problem to improve the accuracy of the solution. This algorithm can be stated as 
follows. 

Solve 

‘[.^]t{Au}« =t (2.111) 

where 

( 2 . 112 ) 


=, {AF} 


(2.113) 



till the convergence criterion 




< tolc 


(2.114) 


is satisfied. The condensed versions of the vectors {f}and^'^^^{R} are obtained 

from their full versions *"‘‘^4/} by the expressions similar to equations 

(2.87) and (2.89). The vector t{R} is called as the unbalanced force vector. The right 
superscript on denotes the configuration on which the integration is to be performed 

to find the external force vector. 

As stated in section 2.11.1, the equation of the first Newton-Raphson iteration (i.e. the 
equation corresponding to =t {AF} ) is solved iteratively using the friction algorithm 

to determine the status (sticking/slipping) of the interface nodes. Corresponding to this 
status, further Newton-Raphson iterations are carried out to reduce the magnitude of the 
unbalanced force vector to zero so as to establish the equilibrium between the external and 
internal forces. 


2.11.3 Arc Length Method 

Since the forging problem is a displacement control problem, there is no known external force. 
Therefore, the denominator in the convergence criterion (2.114) doesn’t exist. In that case, 
one can try to achieve the convergence by making the unbalanced force vector small 

in an absolute sense. But this slows down the rate of convergence to a considerable extent. 
To accelerate the rate of convergence, the arc length method is used in conjunction with 
the modified Newton-Raphson method. In this approach, the unknown nodal reactions at 
the interface are expressed as a linear combination of a set of known nodal forces. Then 
the problem is solved as a force control problem using the set of known nodal forces. The 
unknown coefficients in the linear combination are found from the condition that the vertical 
component of incremental displacement vector at the interface has a specified value. Since 
this solution may not satisfy the nodal equilibrium of internal and external forces, iterations 
need to be performed which are carried out by the modified Newton-Raphson method. Since 
both the nodal reactions as well as the unbalanced force vector change during an iteration, 
the incremental equations contain an additional term compared to eq. (2.111). 

Originally, the arc length method was proposed only for the case of proportional loading 
[83,84]. However, in forging problem, the loading is non-proportional. Therefore, the method 
needs to be modified appropriately to take care of non-proportional loading. The modified 
method is described in this section. 

In iteration, the incremental eq. (2.110) is expressed as^ 


^In this section, the superscripts/subscripts denoting time have been omitted for the sake of convenience. 



(2.115) 


[^]{Au}W = {AF}W + {H}(*-i) 


where 


is the unbalanced force vector given by 


- {/}(*-!) (2.116) 

and the incremental force vector {AF}W is expressed as 

m 

{AF}M = ^(AA)i {P}i (2.117) 

Jk=l 

Here, m is the number of nodes at the interface and {P}fc is a known nodal force vector 
corresponding to a unit (vertically downward) point force at node k. The vectors {P}*, 
k = l,m are called as the basic load vectors. The coefficients (AA)[*^ are unknowns. 

Since equation (2.115) is a linear equation, one can decompose the solution as 

m 

{An}(‘) = ^(AA)UAn(*)^} + (2-118) 

A:=l 

where {Au^*^^}and are obtained as solutions of the following problems: 

[^]{Au<‘)'} = {Pk fc = l,m, (2.119) 

[^]{Au"}<‘> = (2.120) 

Thus, {Au^*'^^}represents the iterative displacement vector due to the basic load vector 
at node k and represents the iterative displacement vector due to the unbalanced 

force vector ' 

To find the complete solution of (2.115), one needs to determine (AA)^*\ k = l,m. They i 
are determined from the following condition. At the nodes on the interface, the vertical or i 
^:-component of the incremental displacement vector is known. Thus, for nodes I = l,m: [ 

m [ 

^(AA)J;^ + Au[f ^ = specified value (2.121) ; 

A;=l : 

These are m equations in m unknowns: (AA)J?\/c = l,m. By solving these equations, 
the unknown (AA)^*^ are calculated. Then the iterative displacements are determined from 
eq (2.118). These displacements are used to find the internal force vector {/}^'^. Finally, the 
external forpe vector at the end (of i^) iteration is found from 

m 

{P}B = 

A;=l 


( 2 . 122 ) 


This force vector is used to find the unbalanced force vector and then to check the con- 
vergence using oq. (2.114). 

While solving eq (2.119,2.120), the upper triangularization of the coefficient matrix [K] 
needs to be done only once. Then, one can use the resolving facility first to perform the Gauss 
elimination operations on different right side vectors and then to find {Aw('')^}and{Aw^^}(") 
by back substitution. 

Equation (2.121) takes a different form for the first iteration than for other iterations. 
These forms are as follows: 

First iteration: 

for the first iteration, = {0} implies that {Au^^}W = {0}. Thus, eq (2.121) 

becomes 

m 

^(AA)^*^Ari[*^^ = specified value (2.123) 

k=l 

Subsequent iterations: 

Since the value of specified incremental displacement does not change during an iteration 
cycle, in subsequent iterations, the right side of eq. (2.121) becomes zero. Thus, 

m 

53(AA)fAai‘>' + Ati{J" = 0 < > 2 (2.124) 

A=1 


2.11.4 Numerical Integration Scheme 


Exact evaluation of integrals appearing in element coefficient matrices and right side vectors 
is not always possible because of the complexity of the integrands. In such cases, it is natural 
to seek numerical evaluation of these integral expressions. Numerical integration involves 
approximation of the integrand by a polynomial of appropriate degree, because the integral 
of a polynomial can be evaluated exactly. 

The most commonly used numerical integration method is Gauss-Legendre integration 
method. The expression for the integral over a master element is given by 



/ +! /*+! 

^ y F((,n)d(,drj 


m n 

-EE^K*- p.l25) 

1=1 j=l 

2.11.5 Divergence Handling Procedures 

The modified Newton-Raphson method diverges in some cases. The following simple but fairly 
effective techniques, for overcoming the divergence are incorporated in the present work. 


1. Under-relaxation: The iterative displacement is scaled by a factor q;„ (0^a„<l). Thus, 
the iterative displacement vector is given by 

(2.126) 

where { Au}^d jg the solution of equation (2.111) and is determined so that the numerical 
method does not diverge. 

2. Line search: In this technique, the displacement vector is updated according to the 
relation [78] 


i+At{-}(i) ^t+At (2.127) 

where ai is determined to minimize the unbalanced force vector. 

A full line search is computationally expensive. Therefore a cheaper technique is to evalu- 
ate unbalance at n different aj values and choose that a; which corresponds to the minimum 
unbalance. This offers a good compromise between computational effort spent in line search 
and effort saved due to convergence in reduced number of iterations. 

3. Incremental displacement cutting: In case the above two procedures fail, incremental 
displacement cutting is performed. 



Chapter 3 

Results and Discussion 


The finite element model of axisymmetric forging process developed in the previous chapter 
has been applied to a number of cases involving various sets of input variables. 

In this chapter, the study is performed on a cylindrical disc. Only upper half of the 
disc is considered for analysis due to symmetry of the geometry and boundary conditions. 
The convergence study is carried out for Dhar’s [46,75] damage variable D and stresses with 
mesh sizes 3x3, 6x6 and 12x12. Two increment sizes 0.1 mm and 0.05 mm are used in the 
convergence study. In section 3.1, some typical results are presented where growth of D 
is discussed along with the hydrostatic stress, circumferential stress ag and axial stress 
distributions. In the last section, the parametric study is carried out to show the effects of 
various input variables. A comparison of bulged profiles is also presented. 

3.1 Typical Results 

In this section, the distribution of Dhar’s damage variable D is closely studied for a typical 
set of input variables. The material properties, friction coefficient and blank size used in this 
study are as follows: 

Material Properties 
Material: Steel SAE 1090, 

Modulus of Elasticity E = 210 GPa, 

Poisson’s ratio p = 0.3, 

Initial yield stre^ {<^y)o — 464.00 MPa, 

Hardening coefficient K =1115.00, 

Hardening Exponent n = 0.19; 

Friction Coefficient; 
fx = 0.05; 


Blank Size: 

Height {H) = 10 nmi. 

Diameter { 0 ) = 10 mm; 

Constants for Dhar’s criterion [46,75]: 

Cd = 1-898 X 10“^ 

m = 9.8 X lO-** MPa-> , 

fl2 = 1.84 MPa->. 

Critical damage value Dc = 0.05. 

A 12x12 elements nonuniform mesh is used and step size taken is 0.05 mm. Three different 
reductions are considered, when the damage variable D reaches the critical value of 0.05 at 
(I) the die-workpiece interface near the free surface of the cylindrical disc, (II) the geometric 
centre of the disc and (III) the middle of the free surface of disc. These reductions are 25%, 
29% and 35% respectively. 

Figures 3.1, 3.2 and 3.3 show the damage distribution at these reductions. Analysis of these 
figures shows that D becomes critical, first at the edge (i.e. at the die-workpiece interface 
near the free surface), then at the centre of the disc. With further deformation, D grows 
radially outwards and then reaches the critical value at the meridian surface (the middle of 
the free surface). By this reduction, D is above 0.05 everywhere in the domain except the 
centre of flat surface of the disc. Equivalent strain distribution at 35% reduction is shown 
in figure 3.4. From this, it is clear that the maximum damage at the edge is due to severe 
deformation and hence high strains here. On the other hand, the strain values are minimum 
at the centre of the flat surface. Thus, micro-cracks will first occur at the edge of the disc, 
then at the centre and finally at the meridian surface. But experimental results [56,57,72,73] 
show that the fractures are seen at the meridian surface. To understand this, hydrostatic 
stress, circumferential stress gb and axial stress distributions are also discussed below. 

Figures 3.5, 3.6 and 3.7 show the distributions of hydrostatic stress at diflferent reductions 
mentioned above. FVom these figures, it is clear that the hydrostatic pressure (negative 
hydrostatic stress) is minimum at the meridian surface of the disc, while it is quite higher at 
the other two locations. Thus, the possibility of a fracture is higher at the meridian surface. 
As the reduction is increased, the hydrostatic pressure becomes less at the meridian surface, 
while it increases at the other two locations suppressing the growth of micro-cracks developed 
there. 

Study of the figures 3.8, 3.9 and 3.10 showing the oe distributions at diflferent reductions 
makes the location of fracture (not micro-cracks) initiation site more clear. It can be seen that 
ae is tensile at the meridian surface, while it is compressive at the top and at the centre of the 
disc. With higher reductions, the value of ae becomes more tensile at the meridian surface, 
while it becomes more compressive at the top and the core of the disc. Thus, micro-cracks 



initiated at the meri<iian surface can grow fast to a fracture. 

At these three reductions, distributions of are also shown (figures 3.11-3.13). The axial 
stre.ss. (T; i.s coinpressivf' everywhere in the domain. It is the minimum compressive at the 
meridian .surface and thi.s shows the maximum chance of crack growth here. Unlike ag, 
continues to hec{)me more compressive everywhere in the domain with increase in reduction. 

Figure 3.14 show.s bulge profiles at the three reductions presented above. The bulging is 
more at higher reduction. 

3.2 Parametric Studies 

In this, section effects of the coefficient of friction and height to diameter ratio on micro-cracks 
occurrence are studied. The nonuniform mesh with 144 elements (12x12) and the increment 
size of 0.05 mm are used, which are same as used in the previous section. The material is also 
the same (SAE 1090 steel). The size of the blank is such that its volume is equal to that of 
the blank of height (ff) and diameter (0) equal to 10 mm each which has been used earlier. 

3.2.1 Effect of Coefficient of Friction 

The analysis is carried out for the coefficient of friction \x — 0.1 and then results are compared 
with that for \x = 0.05 described earlier. The analysis is carried out for H/(i> ratio equal to 1 
[H = ~IQ mm) and for three distinguished reductions as mentioned earlier. 

These reductions in this case are found to be 20%, 25% and 35% when damage D becomes 
critical at the edge, at the centre and at the meridian surface of the disc respectively. A worth 
mentioning point in this case is that D becomes critical at 20% reduction only, while with 
/Li = 0.05, this was at 25% reduction. Thus, with higher friction, micro-cracks initiate at less 
reduction along the die-workpiece interface near the free surface. 

Figures 3.15-3.17 show the growth of D with increasing reduction and the trend is same 
as discussed in the previous section (3.1). Thus micro-cracks first initiate at the edge of disc, 
then at the centre and finally at the meridian surface. At these three reductions, hydrostatic 
stress, erg and cr* distributions are also presented in figures 3.18-3.26. Again the pattern is 
the same as discussed in section 3.1. Hydrostatic stress and cr* are the least compressive and 
ag is the maximum tensile at the meridian surface. At higher reduction, hydrostatic pressure 
becomes more everywhere except the meridian surface where it turns less. The circumferential 
stress ag becomes more tensile and the axial stress less compressive at this location. Here, 
therefore, micro-cracks will grow faster and fracture will occur first. 

Figures 3.3 and 3.17 can be compared for damage at different friction levels. It can be 
noted that the damage D increases everywhere in the domain except the region near the 
centre of the flat face of the disc when the friction is increased. With n = 0.1, becomes 
more compressive everywhere except the region near the meridian surface where it turns less 


coniprossiv(': whih* og hcronu's more compressive at the edge and at the core of the disc and, 
more tensiie at the meridian surface. At higher friction, the hydrostatic pressure increases 
everywhere except the meridian surface where it reduces and becomes nearly zero. This shows 
that the fracture ot'curs at less reduction with higher friction. 

Figure 3.27 shows the bulge profiles for /x = 0.05 and ^ = 0.1, at 35% reduction. The 
bulge is severe with higher friction due to higher resistance to material movement along the 
die-workpiece interface. 

3.2.2 Effect of Height to Diameter Ratio (H/(f> Ratio) 

The effect of height to diameter ratio ratio) is studied by carrying out the analysis for 
the following three cases: (I) H/<f> =0.5, (II) H/(j> =1, (III) H/(t> =2. While varying the H/cj) 
ratio, the volume of the blank is kept constant. The analysis is carried out for the coefficient 
of friction /i =0.05. Three typical reductions are considered when the damage D reaches the 
critical value at three locations in the disc: the edge, centre and meridian surface respectively. 

These reductions are found to be 27.2%, 28.8% and 31.9% for H/4> =0.5; 25.0%, 29.0% 
and 35.0% for H/4> = l and 27.3%, 30.4% and 32.9% for =2. This shows that the micro- 
cracks initiate at less reduction in case of H/(f> =1. Damage D, hydrostatic stress, erg and cr^ 
distributions for H/<t> =0.5 are shown in figures 3.28-3.39, and for H/ej) =2, in figures 3.40- 
3.51. These distributions for H/^ =1 are already shown in figures 3. 1-3.3 and 3.5-3.13. For 
a particular H/4) ratio, the pattern of these distributions and their variation when reduction 
is increased is the same as discussed in section 3.1. 

The cases for various H/<p ratios when D is 0.05 near the origin (disc-centre) are considered 
next. Figures 3.32, 3.6 and 3.44 show that the hydrostatic pressure, ag and at the centre 
of the disc become less compressive with increasing H/4> ratio. Thus, it can be said that the 
possibility of formation of a central cavity is more with the higher Hf ^ ratio. 

Comparison of figures 3.33, 3.7 and 3.45 (when the damage D has reached 0.05 at the 
meridian surface) show that in case of =1, the hydrostatic pressure at the meridian 
surface is less than that in case of =0.5 or HJ<j) =2. Similarly, figures 3.36, 3.10 and 3.48 
show that ag is the maximum tensile for =1. Thus meridian surface fracture will appear 
at less reduction for jfif/0 =1- This is in agreement with the work of Clift et al. [65]. 
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Figure 3.5: Hydrostatic Stress Distribution (in MPa). (Reduction=25.0%, fi=0.05, H/(j)=l, 
H=10.0 mm, (j>=10.0 mm, Mesh Si2e=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.6: Hydrostatic Stress Distribution (in MPa). (Reduction=29.0%, |r=0.05, H/(j)-l, 
H=10.0 mm, 6 = 10.0 mm, Mesh Size?=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.7: Hydrostatic Stress Distribution (in MPa). (Reduction=35.0%, ^=0.05, H/(j)=l, 
H=10.0 mm, <i>=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.9: 00 Distribution (in MPa). (Reduction=29.0%, }x=0.05, H/(j)=l, H=10.0 mm, 
<j)=!0.0 mm, Mesh Size=J2xl2, Increment size=0.05 mm, AISI-1090 steel) 



Figure 3.10: 09 Distribution (in MPa). (Reduction— 35.0%, fi— 0.05, H/(t) 1,H 10.0 mm, 

({>=10.0 mm, Me^ Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 




Figure 3.11; o* Distribution (in MPa). (Reduction=25.0%, n=0.05, H/<j)=l, H=10.0 mm, 
(|>=1 0.0 mm, Mesh Si2e=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.12: (Tz Distribution (in MPa). (Reduction=29.0%, M-05, H/<|)=1, H=10.0 mm, 
(j>=10.0 mm, Mesh Siz©=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.1 3: Distribution (in MPa). (Reduction=35.0%, |a.=0.05, H/(|)=l, H=10.0 mm, 

(j)=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3 14: Bulge Profiles at Different Reductions ()i-0. 05, H/(j) 1,H 10.0 mm, <|) 10.0 mm 
Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.19: Hydrostatic Stress Distribution (in MPa). (Reduction=25.0%, |J.=0.1, H/(j)=l, 
H=i0.0 mm, <|)=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 


Figure 3.20: Hydrostatic Stress Distribution (in MPa). (Reduction=35.0%, n=0.1, H/(l)=l 
H=10.0 mm, f=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.21: ere Distribution (in MPa). (Reduction=20.0%, ^i=0. 1 , H/(j)=l , H=1 0.0 mm, 
(j)=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.22: cq Distribution (in MPa). (Reduction=25.0%, ji=0.1, H/cj)-!, H-10.0 mm, 
<t>=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



439 (max) 


Figure 3.23: ae Distribution (in MPa). (Reduction=35.0%, ^=0.1, H/(f)=l, H=10.0 mm, 
<|>=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.24: Gz Distribution (in MPa). (Reduction=20.0%, n=0.1, H/<j)=l, H=10.0 mm, 
d>«10.0 mm, M^h Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.25: Gz Distribution (in MPa). (Reduction=25.0%, )i=0.1, H/(j)=l, H=10.0 mm, 
(|)=10.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI-1090 steel) 



Figure 3.26: O 2 Distribution (in MPa). (Reduction=35.0%, n=0.1, H/(t)=l, H=10.0 mm, 
6 =\ 0.0 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 




z-axis 



Figure 3.27: Bulge Profiles at Different Friction Levels. 
(Reduction=35.0%, H/<t)=l, H=10.0 mm, (t)=10.0 mm, Mesh Size=12xl2, 
Increment size=0.05 mm, AISI -1090 steel) 
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0.05 (max) 


Figure 3.28: Damage ‘D’ Distribution. (Reduction=27.2%, |x=0.05, H/<i)=0.5, H=6.3 mm, 
<t)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI-1090 steel) 



0.08 (max) 


Figure 3.29: Damage ‘D’ Distribution. (Reduction=28.8%, ^=0.05, H/(j)-0.5, H-6.3 mm, 
(b=12.6 mm. Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.30: Damage ‘D’ Distribution. (Reduction=31.9%, [a,=0.05, H/(j)=0.5, H=6.3 mm, 
<i)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.3 1. : Hydrostatic Stress Distribution (in MPa). (Reduction=27 .2%, ^=0.05, H/(j)-0.5, 
H=6.3 mm, <j>=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.32: Hydrostatic Stress Distribution (in MPa). (Reduction=28.8%, |j,=0.05, H/(j)=0.5, 
H=6.3 mm, <t)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.33: Hydrostatic Stress Distribution (in MPa). (Reduction=31.9%, |x=0.05, H/(|)=0.5, 
H==6.3 mm, ^=\ 2.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.34: 0 e Distribution (in MPa). (Reduction=27.2%, |j,=0.05, H/(j)=0.5, H=6.3 mm, 
(j>=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.35: 0e Distribution (in MPa). (Reduction— 28.8%, |j,— 0.05, Hy(|)— 0.5, H 6.3 mm, 
<j)=12.6 mm. Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.36: oe Distribution (in MPa). (Reduction=3 1 .9%, )a=0.05, H/(j)=0.5, H-6.3 mm, 
(j)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.37: CTz Distribution (in MPa). (Reduction=27.2%, fj.=0.05, H/(t)-0.5, H-6.3 mm, 
^=12.6 mm. Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.38: Distribution (in MPa). (Reduction=28.8%, )j.=0.05, H/(t)=0.5, H=6.3 nun, 

<j)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3 .39: Distribution (in MPa). (Reduction=3 1 .9%, )j.=0.05, H/cj)— 0.5, H-6.3 mm, 

(j)=12.6 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 
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Figure 3.40: Damage ‘D’ Distribution. (Reduction=27.3%, |a,=0.05, H/(t)=2, H=15.8 mm, 
(j)=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.41; Damage ‘D’ Distribution. (Reduction=30.4%, ^i=0.05, H/(j)-2, H-15.8 mm 
*=7,9 mm. Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.42: Damage ‘D’ Distribution. (Reduction=32.9%, )j,=0.05, H/(|)=2, H=15.8 mm, 
4=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.43: Hydrostatic Stress Distribution (in MPa). (Reduction=27.3%, ^=0.05, H/(t)-2, 
H=15.8 mm, 4=7.9 mm. Mesh Size^l2xl2, Increment size=0.05 nam, AISI -1090 steel) 
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Figure 3.44: Hydrostatic Stress Distribution (in MPa). (Reduction=30.4%, p,=0.05, H/(t)-2, 
H=15.8 mm, (j)=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Fieure 3 45’ Hydrostatic Stress Distribution (in MPa). (Reduction 32.9%, p.-0.05, H/(t> 2, 
8 mm,“=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mn, AISI -1090 steel) 




Figure 3.46: ag Distribution (in MPa). (Reduction=27.3%, ji=0.05, H/(t)-2, H-15.8 mm, 
<j)-7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.47: <Je Distribution (in MPa). (Reduction-30.4%, }j,-0.05, H/(|) 2, H 15.8 mm, 
^=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 




Figure 3.48: <jq Distribution (in MPa). (Reduction=32.9%, ^=0.05, H/(t)=2, H=15.8 mm, 
<t)=7.9 mm. Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 



Figure 3.49: az Distribution (in MPa). (Re(iuction=27.3%, ^1=0.05, H/(j)=2, H=15.8 mm, 
<t>=7.9 mm, Mesh Size=12xl2, Incremrat siz^.05 mm, AISI -1090 steel) 





Figure 3.50: Oz Distribution (in MPa). (Reduction=30.4%, |j,=0.05, H/(j)=2, H=15.8 mm, 
<j>=7.9 mm, Mesh Size=12xl2, Increment size=0.05 mm, AISI -1090 steel) 


Figure 3.5 1 : dz Distribution (in MPa). (Reduction-32.9%, p,-0.05, H/<|) 2, H 15.8 mm, 
(J)=7.9 mm, Mesh Size=12xl2, Increment size^.05 mm, AISI -1090 steel) 
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Chapter 4 


Conclusions and Scope for Future 
Work 


The objective of this work is to develop an axisymmetric large deformation elasto-plastic finite 
element code for determination of damage in cold forging process. Coulomb friction model 
is employed at the die-workpiece interface. The Updated Lagrangian formulation, which is 
convenient for handling material and geometric non-linearities, is used. New incremental 
objective stress measure and logarithmic strain measure are used which allow the use of large 
incremental displacement. Forging process is a displacement-control problem. Therefore, 
to accelerate the convergence of the iterative solution scheme, arc length method is used in 
conjunction with the modified Newton Raphson iterative technique. The material is assumed 
to be elastic-plastic strain hardening and yielding according to the von Mises criterion. The 
strain hardening behaviour is modelled by a power law. The inertial and body forces are 
neglected. 

Ductile fracture is often a limiting factor in metal forming processes. It is widely known 
that the ductile fracture is a void nucleation, growth and coalescence phenomenon. A con- 
tinuum damage mechanics model given by Dhar et al. [46,75] is used for the prediction of 
defects. For typical values of input variables, the damage, hydrostatic stress, circumferential 
stress and axial stress distributions are studied. A detailed parametric study of damage is 
carried out to show the effects of the coefficient of friction at the die-work interface and the 
height to diameter ratio of the workpiece. 

The following conclusions can be drawn from the present study: 

• The micro-cracks first initiate along the die-workpiece interface region near the free 
surface, then at the centre of the workpiece and subsequently at the meridian surface. 
The amount of damage continues to grow with reduction. 

• The fracture will appear first at the meridian surface since the hydrostatic stress and 
axial stress are the minimum compressive and the circumferential stress is the maximum 
tensile at this location. Thus, micro-cracks here grow faster to a fracture. 
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• .-M mgiu-i uR iiuii, imcro-cracKs inmate at less reduction, i he meridian suriace tracture 
will also ocTur early with higher friction since the circumferential stress is more tensile 
anci hydrostatic stress, less compressive at this location. Bulging is more severe when 
the friction at die- work interface is higher. 

• As compared with the case of height to diameter ratio equal to one, blanks with height to 
diameter ratio le.ss than or greater than unity can sustain more reductions without micro- 
cracks and fracture. With higher height to diameter ratio, the chances of formation of 
a central cavity are more because of less compressive hydrostatic stress. 

TIh’ following points may be considered as an extension of the present study: 

• Dependence of material behaviour on strain rate (visco-plasticity) and temperature. 

• Analysis and prediction of defects in forging process during unloading and due to residual 
stresses. 

• Determination of the properties for other materials to fit the criterion of void nucleation, 
growth and coalescence. 

• Preprocessor and postprocessor for generation of input data and graphical representation 
of results. 
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